eV=[0.64 0.77 0.89 1.02 1.14 1.26 1.39 1.51 1.64 1.76 1.88 2.01 2.13 2.26 2.38 2.50 2.63 2.75 2.88 3.00 3.12 3.25 3.37 3.50 3.62 3.74 3.87 3.99 4.12 4.24 4.36 4.49 4.61 4.74 4.86 4.98 5.11 5.23 5.36 5.48 5.60 5.73 5.85 5.98 6.10 6.22 6.35 6.47 6.60]';
% gold
n=[0.92 0.56 0.43 0.35 0.27 0.22 0.17 0.16 0.14 0.13 0.14 0.21 0.29 0.43 0.62 1.04 1.31 1.38 1.45 1.46 1.47 1.46 1.48 1.50 1.48 1.48 1.54 1.53 1.53 1.49 1.47 1.43 1.38 1.35 1.33 1.33 1.32 1.32 1.30 1.31 1.30 1.30 1.30 1.30 1.33 1.33 1.34 1.32 1.28]';
k=[13.78 11.21 9.519 8.145 7.150 6.350 5.663 5.083 4.542 4.103 3.697 3.272 2.863 2.455 2.081 1.833 1.849 1.914 1.948 1.958 1.952 1.933 1.895 1.866 1.871 1.883 1.898 1.893 1.889 1.878 1.869 1.847 1.803 1.749 1.688 1.631 1.577 1.536 1.497 1.460 1.427 1.387 1.350 1.304 1.277 1.251 1.226 1.203 1.188]';
e1=n.^2-k.^2;
e2=2*n.*k;
% silver
% n=[0.24 0.15 0.13 0.09 0.04 0.04 0.04 0.04 0.03 0.04 0.05 0.06 0.05 0.06 0.05 0.05 0.05 0.04 0.04 0.05 0.05 0.05 0.07 0.10 0.14 0.17 0.81 1.13 1.34 1.89 1.41 1.41 1.38 1.35 1.38 1.31 1.30 1.28 1.28 1.26 1.25 1.22 1.20 1.18 1.15 1.14 1.12 1.10 1.07]';
% k=[14.08 11.85 10.10 8.828 7.795 6.992 6.312 5.727 5.242 4.838 4.483 4.152 3.858 3.586 3.324 3.093 2.869 2.657 2.462 2.275 2.070 1.864 1.657 1.419 1.142 0.829 0.392 0.616 0.964 1.161 1.264 1.331 1.372 1.387 1.393 1.389 1.378 1.367 1.357 1.344 1.342 1.336 1.325 1.312 1.296 1.277 1.255 1.232 1.212]';


%lambda_nm=100:2:1500;
%lambda=1e-9*lambda_nm;
hbar=jcb_hbar;
c=jcb_c;
q=1.602176487e-19;

w=eV2w(0.5:0.001:6);

wp=13.8e15;
epsinf=8.7;
tau=1.0e-14;
gamma=1/tau;
epsilon_d=epsinf-(wp^2)./(w.^2+1i*gamma.*w);

n_d=real(sqrt(epsilon_d));
k_d=imag(sqrt(epsilon_d));

%%
figure(1);
clf
ax(1)=subplot(1,2,1);
h1=plot(ax(1),w2eV(w),real(epsilon_d));
ax(2)=subplot(1,2,2);
h2=plot(ax(2),w2eV(w),imag(epsilon_d));
set(h1,'color','b')
set(h2,'color','r')
%set(ax,'YColor','k')
set(get(ax(1),'YLabel'),'String','\Ree[\epsilon(\omega)]');
set(get(ax(2),'YLabel'),'String','\Imm[\epsilon(\omega)]');
set([get(ax(1),'XLabel') get(ax(2),'XLabel')],'String','\lambda/nm');
set([ax get(ax(1),'Xlabel') get(ax(2),'Xlabel') get(ax(1),'Ylabel') ... 
    get(ax(2),'Ylabel') get(ax(1),'title') get(ax(2),'title')],...
    'Fontsize',8,'FontName','Times New Roman')

set(ax(1),'Ylim',[-30 5])
set(ax(2),'Ylim',[0 7])
set(ax,'xlim',[0.5 6],'ytickmode','auto')
set(get(ax(1),'title'),'string','\Ree[\epsilon(\omega)]')
set(get(ax(2),'title'),'string','\Imm[\epsilon(\omega)]')
%legend({'\Ree[\epsilon(\omega)]','\Imm[\epsilon(\omega)]'})


%set(ax(2),'xtick',[])
xtick=w2eV(w2lam([2000 1000 700 500 400 300 250 210]*1e-9));
set(ax(1),'xticklabel',num2str(round(1e9*w2lam(eV2w(xtick')))),...
          'xtick',xtick)
set(ax(2),'xticklabel',num2str(round(1e9*w2lam(eV2w(xtick')))),...
          'xtick',xtick)

% add J&C data points
hold(ax(1),'on')
hold(ax(2),'on')
plot(ax(1),eV,e1,'b.','linestyle','none');
plot(ax(2),eV,e2,'r.','linestyle','none');
plot([0.5 6],[0 0],'k')
hold(ax(1),'on')
hold(ax(2),'on')

set(ax(1),'position',[0.07 0.15 0.42 0.75]);
set(ax(2),'position',[0.57 0.15 0.42 0.75]);

% print figure
set(gcf,'color','white');
set(gcf,'PaperPositionMode','auto','units','centimeters',...
       'PaperSize',[16 6],'PaperPosition',[0 0 16 6]);
print('-dpdf','-r1000','DrudeDielectric');
print('-dpng','-r1000','DrudeDielectric');

%% model with lorentz terms

wp=13.8e15;
epsinf=11;
tau=1.2e-14;
gamma=1/tau;
epsilon_d=epsinf-(wp^2)./(w.^2+1i*gamma.*w);

threshold=2.8;
Aj =  [1 1 0.8*[1 0.8 1] 0.8*[1 1 1 1 1] 0.5*[1 1 1 1 1 1 1 1 1] ]*0.4*(45e14)^2;
wj = eV2w([2.8000 3.0000 3.2000 3.4000 3.6000 3.8000 4.0000 4.2000 4.4000...
           4.6000 4.8000 5.0000 5.2000 5.4000 5.6000 5.8000 6.0000 6.2000...
           6.4000 6.6000]);%linspace(w2lam(210e-9),w2lam(450e-9),10);

gammaj = 1e15*1*[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ];%[1 1 1 1 1 1 1 1 1 1]*1e15;
eps_l=cell(size(Aj));
hold on
y=delta_eps_inf+epsilon_d;

for ii=1:length(Aj)
    eps_l{ii} = Aj(ii)./((wj(ii)^2 - w.^2)-1i*gammaj(ii)*w);
    y=y+eps_l{ii};
end

set(h1,'xdata',w2eV(w),'ydata',real(y));
set(h2,'xdata',w2eV(w),'ydata',imag(y))

% print figure
print('-dpdf','-r1000','DrudeLorentzDielectric');
print('-dpng','-r1000','DrudeLorentzDielectric');

%%
% figure(1);
% clf
% [ax h1 h2]=plotyy(w2eV(w),n_d,w2eV(w),k_d);
% set(h1,'color','b')
% set(h2,'color','r')
% set(ax,'YColor','k')
% set(get(ax(1),'YLabel'),'String','n(\omega)')
% set(get(ax(2),'YLabel'),'String','\kappa(\omega)')
% set([ax get(ax(1),'Xlabel') get(ax(1),'Ylabel') get(ax(2),'Ylabel')],...
%     'Fontsize',12,'FontName','Times New Roman')
% set(ax(1),'Ylim',[0 2])
% set(ax(2),'Ylim',[0 20])
% set(ax,'ytickmode','auto')
% set(ax,'xlim',[0.5 4],'box','off','ytickmode','auto')
% title('n,\kappa')
% xlabel('2\pih\omega /eV')
% set(ax(2),'xtick',[])
% legend({'n(\omega)','\kappa(\omega)'})
% 
% % add J&C data points
% hold(ax(1),'on')
% hold(ax(2),'on')
% plot(ax(1),eV,n,'bx','linestyle','none')
% plot(ax(2),eV,k,'r+','linestyle','none')
% hold(ax(1),'on')
% hold(ax(2),'on')

%%
% set(gcf,'color','white');
% set(gcf,'PaperPositionMode','auto','units','pixels','position',[50 50 xsize-1 xsize*3/4]);
% set([waxis laxis baxis],'units','normalized');
% set([waxis laxis baxis],'position',[0.07 0.075 0.9 0.85]);
% print('-dpdf','-r0',picfile);
